/**
 * @file Element.h
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2020-11-08
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _ELEMENT_H_
#define _ELEMENT_H_
#include "Dofs.h"
#include "GaussianPoint.h"
#include <assert.h>
#define TEMPLATE template<Dimension DIM>
/**
 * @brief 有限元计算过程中的基类单元，用于派生出各种不同单元；
 *        Element是TemplateElement的缩写，表示模板单元或者成为referernce element;
 */
TEMPLATE
class Element
{
typedef std::vector<Dofs<DIM> > DofsList;
typedef int DofsNum;
public:
    /**
     * @brief Construct a new Element object
     * 
     */
    Element();
    /**
     * @brief Destroy the Element object
     * 
     */
    ~Element();
    /**
     * @brief Set the Dofs List object by std::vector<Dofs<DIM> >
     * 
     */
    virtual void SetDofsList(DofsList) = 0 ; 
    /**
     * @brief Get the Dofs List object
     * 
     * @param i 指定序号；
     * @return Dofs<DIM>&　自由度的引用 
     */
    virtual Dofs<DIM>& GetDofsList(int i) = 0 ;
    /**
     * @brief 返回本单元里第ｉ个自由度对应的全局编号；
     * 
     * @param i 指定序号；
     * @return int 全局编号；
     */
    virtual int NdIdx(int i)  = 0 ;
    /**
     * @brief 返回本单元对应第i个高斯点坐标信息；
     * 
     * @param i 指定序号
     * @return Real* 存储高斯点的数组
     */
    virtual Real* GaussionPoint(int i) = 0 ;
    /**
     * @brief 返回本单元对应第i个高斯点权重信息；
     * 
     * @param i 指定序号
     * @return Real 高斯点对应权重的数值；
     */
    virtual Real GaussionWeight(int i) = 0 ;
    /**
     * @brief 该类型单元对应的基函数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i　基函数序数
     * @return Real 该基函数的值；
     */
    virtual Real phi(double xi, double eta, int i)= 0 ;
    /**
     * @brief 该类型单元对应的基函数对xi的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return Real 该偏导数的值
     */
    virtual Real phi_xi(double xi, double eta, int i)  = 0;
    /**
     * @brief 该类型单元对应的基函数对eta的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return Real 该偏导数的值
     */
    virtual Real phi_eta(double xi, double eta, int i)  = 0;
    /**
     * @brief 任意单元到该模板单元仿射变换对应的雅克比矩阵行列式在(xi,eta)处的值；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @return Real 值
     */
    virtual Real det_Jacobi(double xi, double eta) = 0 ;
    /**
     * @brief 基函数对x的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return Real 值
     */
    virtual Real phi_x(double xi, double eta, int i)  = 0 ;
    /**
     * @brief 基函数对y的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return Real 值
     */
    virtual Real phi_y(double xi, double eta, int i)  = 0 ;
    /**
     * @brief 计算并返回第ｉ个基函数在xi，eta点处的梯度向量；
     * 
     * @param xi 
     * @param eta 
     * @param i 
     * @return std::vector<Real> 梯度向量；
     */
    virtual std::vector<Real> gradient(double xi,double eta,int i) = 0;
    /**
     * @brief 标准单元中的点在真实单元的x坐标
     * 
     * @param xi 坐标１
     * @param eta 坐标2
     * @return Real 值
     */
    
    virtual Real Global_x(double xi, double eta)  = 0 ;
    /**
     * @brief 标准单元中的点在真实单元的y坐标
     * 
     * @param xi 坐标１
     * @param eta 坐标2
     * @return Real 值
     */
    virtual Real Global_y(double xi, double eta)  = 0;
    /**
     * @brief 该单元自由度个数；
     * 
     * @return int 
     */
    virtual int n_Dofs() = 0;
    /**
     * @brief 高斯点个数；
     * 
     * @return int 
     */
    virtual int n_GaussPnt() = 0;
    /**
     * @brief 计算并返回该单元的面积或体积；
     * 
     * @return double 
     */
    virtual double volume() = 0;
    protected:
    virtual Real xi_x(double xi, double eta)  = 0;
    virtual Real xi_y(double xi, double eta)  = 0 ;
    virtual Real eta_x(double xi, double eta)  = 0 ;
    virtual Real eta_y(double xi, double eta) = 0 ;
    
};
TEMPLATE
Element<DIM>::Element()
{
}
TEMPLATE
Element<DIM>::~Element()
{
}
#undef TEMPLATE

/***************************************************************/
/****************** INTERVAL 1D ELEMENT ************************/
/***************************************************************/
class Interval_1_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
public:
    Interval_1_Element();
    ~Interval_1_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i); 
    virtual Real* GaussionPoint(int i)
    {
	return _GaussInfo.GaussianPointInfo(i);
    };
    virtual Real GaussionWeight(int i)
    {
	return _GaussInfo.GaussianWeightInfo(i);
    };               
    virtual Real phi(double xi, double eta, int i);
    virtual Real phi_xi(double xi, double eta, int i);
    virtual Real phi_eta(double xi, double eta, int i);
    virtual Real det_Jacobi(double xi, double eta);
    virtual Real xi_x(double xi, double eta);
    virtual Real xi_y(double xi, double eta);
    virtual Real eta_x(double xi, double eta);
    virtual Real eta_y(double xi, double eta);
    virtual Real phi_x(double xi, double eta, int i);
    virtual Real phi_y(double xi, double eta, int i);
    virtual std::vector<Real> gradient(double xi,double eta,int i);
    virtual Real Global_x(double xi, double eta);
    virtual Real Global_y(double xi, double eta) ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual int n_Dofs()
    {
	return 2;
    };
    virtual int n_GaussPnt()
    {
	return _GaussInfo.GaussianPointNumInfo();
    };
private:
    DofsList _dofslist;
    Interval_Element_GaussianInfo _GaussInfo;
};

Interval_1_Element::Interval_1_Element()
{
    _dofslist.resize(2);
    _GaussInfo.LoadGaussianInfo(3);
}

Interval_1_Element::~Interval_1_Element(){};

void Interval_1_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;
}

Dofs<2>& Interval_1_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Interval_1_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 0.5 - 0.5 * xi;
    case 2:
	return 0.5 + 0.5 * xi;
    }
    return 0;
}

double Interval_1_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -0.5;
    case 2:
	return 0.5;
    }
    return 0;
}

double Interval_1_Element::phi_eta(double xi, double eta, int i)
{
    return 0;
}

double Interval_1_Element::det_Jacobi(double xi, double eta)
{
    /*double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= 2; i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);*/
    return (_dofslist[1][0] - _dofslist[0][0])/2.0;
}

double Interval_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Interval_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Interval_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Interval_1_Element::eta_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Interval_1_Element::phi_x(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Interval_1_Element::phi_y(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<Real> Interval_1_Element::gradient(double xi,double eta,int i)
{
    std::vector<Real> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}

double Interval_1_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Interval_1_Element::Global_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 2; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Interval_1_Element::NdIdx(int i)
{
    return _dofslist[i-1].GlobalIndex();
}

/***************************************************************/
/********************** Q1 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_1_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_1_Element();
    ~Quadrilateral_1_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual Real* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual Real GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual Real phi(double xi, double eta, int i);
    virtual Real phi_xi(double xi, double eta, int i)  ;
    virtual Real phi_eta(double xi, double eta, int i)  ;
    virtual Real det_Jacobi(double xi, double eta)   ;
    virtual Real xi_x(double xi, double eta)  ;
    virtual Real xi_y(double xi, double eta)   ;
    virtual Real eta_x(double xi, double eta)   ;
    virtual Real eta_y(double xi, double eta)   ;
    virtual Real phi_x(double xi, double eta, int i)  ;
    virtual Real phi_y(double xi, double eta, int i)  ;
    virtual std::vector<Real> gradient(double xi,double eta,int i) ;
    virtual Real Global_x(double xi, double eta) ;
    virtual Real Global_y(double xi, double eta)  ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual int n_Dofs(){return 4;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};

Quadrilateral_1_Element::Quadrilateral_1_Element()
{
    _dofslist.resize(4);
    _GaussInfo.LoadGaussianInfo(3);
};

Quadrilateral_1_Element::~Quadrilateral_1_Element(){};
void Quadrilateral_1_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;

}
Dofs<2>& Quadrilateral_1_Element::GetDofsList(int i)
{
    return _dofslist[i];
}
double Quadrilateral_1_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (xi-1)*(eta-1)/4.0;
    case 2:
	return -(xi+1)*(eta-1)/4.0;
    case 3:
	return (xi+1)*(eta+1)/4.0;
    case 4:
	return -(xi-1)*(eta+1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_xi(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return (eta-1)/4.0;
    case 2:
	return -(eta-1)/4.0;
    case 3:
	return (eta+1)/4.0;
    case 4:
	return -(eta+1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_eta(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return (xi-1)/4.0;
    case 2:
	return -(xi+1)/4.0;
    case 3:
	return (xi+1)/4.0;
    case 4:
	return -(xi-1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= 4; i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::phi_x(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_1_Element::phi_y(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<Real> Quadrilateral_1_Element::gradient(double xi,double eta,int i)
{
    std::vector<Real> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Quadrilateral_1_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Quadrilateral_1_Element::Global_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Quadrilateral_1_Element::NdIdx(int i)
{
    switch(i)
    {
    case 1:
	return _dofslist[0].GlobalIndex();
    case 2:
	return _dofslist[1].GlobalIndex();
    case 3:
	return _dofslist[2].GlobalIndex();
    case 4:
	return _dofslist[3].GlobalIndex();
    }
    return 0;
}

/***************************************************************/
/********************** P1 ELEMENT *****************************/
/***************************************************************/

class Triangular_1_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Triangular_1_Element();
    ~Triangular_1_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual Real* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual Real GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};                      
    virtual Real phi(double xi, double eta, int i);
    virtual Real phi_xi(double xi, double eta, int i)  ;
    virtual Real phi_eta(double xi, double eta, int i)  ;
    virtual Real det_Jacobi(double xi, double eta)   ;
    virtual Real xi_x(double xi, double eta)  ;
    virtual Real xi_y(double xi, double eta)   ;
    virtual Real eta_x(double xi, double eta)   ;
    virtual Real eta_y(double xi, double eta)   ;
    virtual Real phi_x(double xi, double eta, int i)  ;
    virtual Real phi_y(double xi, double eta, int i)  ;
    virtual std::vector<Real> gradient(double xi,double eta,int i) ;
    virtual double volume(){return det_Jacobi(0,0) / 2.0;};
    virtual Real Global_x(double xi, double eta) ;
    virtual Real Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 3;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_1_Element::Triangular_1_Element()
{
    _dofslist.resize(3);
    _GaussInfo.LoadGaussianInfo(4);
}

Triangular_1_Element::~Triangular_1_Element()
{
};

void Triangular_1_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_1_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_1_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1 - xi - eta);
    case 2:
	return xi;
    case 3:
	return eta;
    }
    return 0;
}

double Triangular_1_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -1;
    case 2:
	return 1;
    case 3:
	return 0;
    }
    return 0;
}

double Triangular_1_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -1;
    case 2:
	return 0;
    case 3:
	return 1;
    }
    return 0;
}

double Triangular_1_Element::det_Jacobi(double xi, double eta)
{
    return (_dofslist[1][0] - _dofslist[0][0])*(_dofslist[2][1] - _dofslist[0][1])
	-(_dofslist[1][1] - _dofslist[0][1])*(_dofslist[2][0]-_dofslist[0][0]);
}

double Triangular_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_1_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<Real> Triangular_1_Element::gradient(double xi,double eta,int i)
{
    std::vector<Real> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_1_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_1_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_1_Element::NdIdx(int i)
{
    switch(i)
    {
    case 1:
	return _dofslist[0].GlobalIndex();
    case 2:
	return _dofslist[1].GlobalIndex();
    case 3:
	return _dofslist[2].GlobalIndex();
    }
    return 0;
}

/***************************************************************/
/********************** Q2 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_2_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_2_Element();
    ~Quadrilateral_2_Element(){};
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual Real* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual Real GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual Real phi(double xi, double eta, int i);
    virtual Real phi_xi(double xi, double eta, int i)  ;
    virtual Real phi_eta(double xi, double eta, int i)  ;
    virtual Real det_Jacobi(double xi, double eta)   ;
    virtual Real xi_x(double xi, double eta)  ;
    virtual Real xi_y(double xi, double eta)   ;
    virtual Real eta_x(double xi, double eta)   ;
    virtual Real eta_y(double xi, double eta)   ;
    virtual Real phi_x(double xi, double eta, int i)  ;
    virtual Real phi_y(double xi, double eta, int i)  ;
    virtual std::vector<Real> gradient(double xi,double eta,int i) ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual Real Global_x(double xi, double eta) ;
    virtual Real Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 9;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};
Quadrilateral_2_Element::Quadrilateral_2_Element()
{
    _dofslist.resize(9);
    _GaussInfo.LoadGaussianInfo(5);
}
void Quadrilateral_2_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;
}

Dofs<2>& Quadrilateral_2_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Quadrilateral_2_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*xi*eta*(1.0-xi)*(1.0-eta);
    case 2:
	return -(1/4.0)*xi*eta*(1.0+xi)*(1.0-eta);
    case 3:
	return (1/4.0)*xi*eta*(1.0+xi)*(1.0+eta);
    case 4:
    return -(1/4.0)*xi*eta*(1.0-xi)*(1.0+eta);
    case 5:
    return -(1/2.0)*(1.0-eta)*eta*(1.0-xi*xi);
    case 6:
    return (1/2.0)*(1.0+xi)*xi*(1.0-eta*eta);
    case 7:
    return (1/2.0)*(1.0+eta)*eta*(1.0-xi*xi);
    case 8:
    return -(1/2.0)*(1.0-xi)*xi*(1.0-eta*eta);
    case 9:
    return (1.0-xi*xi)*(1.0-eta*eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*(1.0-2.0*xi)*eta*(1.0-eta);
    case 2:
	return -(1/4.0)*(1.0+2.0*xi)*eta*(1.0-eta);
    case 3:
	return (1/4.0)*(1.0+2.0*xi)*eta*(1.0+eta);
    case 4:
    return -(1/4.0)*(1.0-2.0*xi)*eta*(1.0+eta);
    case 5:
    return xi * (1.0-eta) * eta;
    case 6:
    return (1/2.0)*(1.0 + 2.0*xi) * (1.0-eta*eta);
    case 7:
    return -xi * (1.0+eta) * eta;
    case 8:
    return -(1/2.0)*(1.0 - 2.0*xi) * (1.0-eta*eta);
    case 9:
    return -2.0 * xi * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*xi*(1.0 - xi)*(1.0-2.0*eta);
    case 2:
	return - (1/4.0)*xi*(1.0 + xi)*(1-2.0*eta);
    case 3:
	return (1/4.0)*xi*(1.0 + xi)*(1.0+2.0*eta);
    case 4:
    return - (1/4.0)*xi*(1.0 - xi)*(1+2.0*eta);
    case 5:
    return - (1/2.0)*(1.0 - xi * xi)*(1.0-2.0*eta);
    case 6:
    return - xi * (1.0+xi) * eta;
    case 7:
    return (1/2.0)*(1.0 - xi * xi)*(1.0+2.0*eta);
    case 8:
    return xi * (1.0-xi) * eta;
    case 9:
    return -2.0 * (1.0 - xi * xi) * eta;
    }
    return 0;
}

double Quadrilateral_2_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}
double Quadrilateral_2_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_2_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<Real> Quadrilateral_2_Element::gradient(double xi,double eta,int i)
{
    std::vector<Real> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}

double Quadrilateral_2_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= 9; i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_2_Element::Global_x(double xi, double eta) 
{
    /*
    double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }*/
    double a = 0;
    a += _dofslist[0][0] * 0.25*(1 - xi)*(1 - eta);
    a += _dofslist[1][0] * 0.25*(1 + xi)*(1 - eta);
    a += _dofslist[2][0] * 0.25*(1 + xi)*(1 + eta);
    a += _dofslist[3][0] * 0.25*(1 - xi)*(1 + eta);
    return a;
}

double Quadrilateral_2_Element::Global_y(double xi, double eta)
{
    /*double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }*/
    double a = 0;
    a += _dofslist[0][1] * 0.25*(1 - xi)*(1 - eta);
    a += _dofslist[1][1] * 0.25*(1 + xi)*(1 - eta);
    a += _dofslist[2][1] * 0.25*(1 + xi)*(1 + eta);
    a += _dofslist[3][1] * 0.25*(1 - xi)*(1 + eta);
    return a;
}

int Quadrilateral_2_Element::NdIdx(int i)
{
    return _dofslist[i - 1].GlobalIndex();
}

/***************************************************************/
/********************** P2 ELEMENT *****************************/
/***************************************************************/

class Triangular_2_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Triangular_2_Element();
    ~Triangular_2_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual Real* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual Real GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};                      
    virtual Real phi(double xi, double eta, int i);
    virtual Real phi_xi(double xi, double eta, int i)  ;
    virtual Real phi_eta(double xi, double eta, int i)  ;
    virtual Real det_Jacobi(double xi, double eta)   ;
    virtual Real xi_x(double xi, double eta)  ;
    virtual Real xi_y(double xi, double eta)   ;
    virtual Real eta_x(double xi, double eta)   ;
    virtual Real eta_y(double xi, double eta)   ;
    virtual Real phi_x(double xi, double eta, int i)  ;
    virtual Real phi_y(double xi, double eta, int i)  ;
    virtual std::vector<Real> gradient(double xi,double eta,int i) ;
    virtual double volume(){return det_Jacobi(0,0);};
    virtual Real Global_x(double xi, double eta) ;
    virtual Real Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 6;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_2_Element::Triangular_2_Element()
{
    _dofslist.resize(6);
    _GaussInfo.LoadGaussianInfo(5);
}

Triangular_2_Element::~Triangular_2_Element()
{
};

void Triangular_2_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_2_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_2_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1 - xi - eta)*(2 *(1 - xi -eta) - 1);
    case 2:
	return xi * (2 * xi - 1);
    case 3:
	return eta * ( 2 * eta - 1);
    case 4:
    return 4 * xi * (1 - xi - eta);
    case 5:
    return 4 * xi * eta;
    case 6:
    return 4 * eta *(1 - xi - eta);
    }
    return 0;
}

double Triangular_2_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 1 - 4*(1- xi - eta);
    case 2:
	return 4 * xi - 1;
    case 3:
	return 0;
    case 4:
    return 4 - 8 * xi - 4 * eta;
    case 5:
    return 4 * eta;
    case 6:
    return - 4 * eta;
    }
    return 0;
}

double Triangular_2_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 1 - 4*(1- xi - eta);
    case 2:
	return 0;
    case 3:
	return 4 * eta - 1;
    case 4:
    return - 4 * xi;
    case 5:
    return  4 * xi;
    case 6:
    return 4 - 4 * xi - 8 * eta;
    }
    return 0;
}

double Triangular_2_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Triangular_2_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_2_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<Real> Triangular_2_Element::gradient(double xi,double eta,int i)
{
    std::vector<Real> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_2_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_2_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_2_Element::NdIdx(int i)
{
    
    return _dofslist[i - 1].GlobalIndex();
}

template<typename vector_type1,typename vector_type2>
Real InnerProuduct(const vector_type1 v1, const vector_type2 v2)
{
    return 0.0;
}
template<>
Real InnerProuduct<std::vector<Real>,std::vector<Real>>(const std::vector<Real> v1, const std::vector<Real> v2)
{
    assert(v1.size() == v2.size());
    double ans = 0.0;
    for(int i = 0;i < v1.size();i++)
        ans += v1[i] * v2[i];
    return ans;
}
#else

#endif
